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An explicit finite-difference scheme is presented for solving the two-dimensional Biot equations of 
' poroelasticity across the full range of frequencies. The key difficulty is to discretize the Johnson- 

Koplik-Dashen (JKD) model which describes the viscous dissipations in the pores. Indeed, the time- 
domain version of Biot-JKD model involves order 1/2 shifted fractional derivatives which amounts to 
£\j ■ a time convolution product. To avoid storing the past values of the solution, a diffusive representation 

of fractional derivatives is used: the convolution kernel is replaced by a finite number of memory 
variables that satisfy local-in-time ordinary differential equations. The coefficients of the diffusive 
representation follow from an optimization procedure of the dispersion relation. Then, various 
methods of scientific computing are applied: the propagative part of the equations is discretized using 
a fourth-order ADER scheme, whereas the diffusive part is solved exactly. An immersed interface 
method is implemented to discretize the geometry on a Cartesian grid, and also to enforce the jump 
conditions at interfaces. Numerical experiments are proposed in various realistic configurations. 
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Two frequency regimes have to be distinguished when dealing with poroelastic waves. In the low-frequency range (LF), 1 
the viscous efforts are proportional to the relative velocity of the fluid and the solid matrix. In the high-frequency range 
(HF), modeling the dissipation is a more delicate task: Biot first presented an expression for particular pore geometries^ In 
1987, Johnson-Koplik-Dashcn (JKD) 3 published a general expression for the dissipation in the case of random pores. The 
\q \ viscous efforts depend in this model on the square root of the frequency. In the time domain, time fractional derivatives 
■ are then introduced, which involves convolution products with singular kernels. Transient analytical solutions of Biot-JKD 
' model have been derived in academic geometries and homogeneous medial When dealing with more complex geometries, 
. numerical methods are required. 
i—\ ' Many time-domain simulation methods have been developed in the LF regime: see Ref. [5| and the introduction to 
(N '■ Ref. H for general reviews. In the HF regime, the fractional derivatives greatly complicate the numerical modeling of the 
Biot-JKD equations. The past values of the solution are indeed required in order to evaluate the convolution products, 
which means that the time evolution of the solution must be stored. It greatly increases the memory requirements and 
makes large-scale simulations impossible. To our knowledge, two approaches have been proposed so far in the literature. 
The first approach consists in discretizing the convolution products^ and the second one is based on the use of a diffusive 
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^ ■ representation of the fractional derivative^ In the latter approach, the convolution product is replaced by a continuum of 
diffusive variables - or memory variables - satisfying local differential equations^ This continuum is then discretized using 
appropriate quadrature formulas, resulting in the Biot-DA (diffusive approximation) model. 

In Ref. [ToL the diffusive approach was followed in a one-dimensional homogeneous configuration. Compared with Ref. @, 
important improvements were introduced: good representation of the viscous dissipation on the whole range of frequencies; 
optimization of the model on the range of interest; estimation of the computational effort in terms of the required accuracy. 
The goal of the present work is to extend the algorithm to two-dimensional heterogeneous porous media with interfaces. 

This paper is organized as follows. The original Biot-JKD model is briefly outlined in section fll] and the principles 
underlying the diffusive representation of fractional derivatives are described. The decrease of energy and the dispersion 
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analysis are addressed. In section IIII1 the method used to discretize the diffusive model is presented: the diffusive 
approximation thus obtained is easily treatable by computers. Following a similar approach than in viscoelasticity, 11 
the coefficients of the model are determined using an optimization procedure in the frequency range of interest, giving 
an optimum number of additional computational arrays. The numerical modeling is addressed in section IIV1 where the 
equations of evolution are split into two parts: a propagative part, which is discretized using a fourth-order scheme for 
hyperbolic equations, and a diffusive part, which is solved exactly. An immersed interface method accounts for the jump 
conditions and for the geometry of the interfaces, preventing from the usual limitations of finite differences on a Cartesian 
grid. Numerical experiments performed with realistic values of the physical parameters are presented in section [V] In 
section I VII a conclusion is drawn and some futures lines of research are given. 



II. PHYSICAL MODELING 
A. Biot model 

The Biot model describes the propagation of mechanical waves in a macroscopic porous medium consisting of a solid 
matrix saturated with a fluid circulating freely through the pores >i i 12 i 13 It is assumed that 

• the wavelengths are large in comparison with the diameter of the pores; 

• the amplitude of the perturbations is small; 

• the elastic and isotropic matrix is completely saturated with a single fluid phase; 

• the thermo-mechanical effects are neglected. 

This model involves 10 physical parameters: the density pf and the dynamic viscosity rj of the fluid; the density p s and 
the shear modulus p of the elastic skeleton; the porosity < <f> < 1, the tortuosity a > 1, the absolute permeability k, 
the Lame coefficient Xf and the two Biot's coefficients f3 and m of the saturated matrix. The following notations are 
introduced 

Pw = -rPf, P = 4>Pf + (1 - 4>)p s , X = PPw~p}>0, \ = \ f -mf3 2 . (1) 



Taking u s and Uf to denote the solid and fluid displacements, the unknowns in 2D are the elastic velocity v s = the 
filtration velocity w = ^j- — 4>-^j (uf — u s ), the elastic stress tensor <x, and the acoustic pressure p. The constitutive 
laws are 

a = (A/ tr e — j3 m £) I + 2 /ie, 

(2) 

p = m (— /3tre + £), 

where e = i (Vu s + Vu s T ) is the strain tensor and £ = — V.W is the rate of fluid change. The symmetry of er implies 
compatibility conditions between spatial derivatives of e, leading to the Beltrami-Michell equation: 14 



d 2 _ a 9 2 a xx d 2 a yy d 2 p d 2 a xx d 2 a yy d 2 p 

-5 — 7T~ — t>0 a 9 r V\ » 7, r "2 tt^t + "I ^ 9 r u ,-, 9 r 02 ttt , 

ax ay ox z ox z ox z oy z oy z oy z 

n ^0 n X a + 2p flfi 



(3) 



4(A + m)' " 4(A + /i)' 2 2(A + /i)' 
On the other hand, the conservation of momentum yields 

<9v s 9w 

dv s dw 77 
Ps -qT- + Pw-gj: + ~ F *™ ' = -vp, 



(4) 



where * is the convolution product in time. The second equation of Q is a generalized Darcy law. The quantity F * w 
denotes the viscous dissipation induced by the relative motion between the fluid and the elastic skeleton. 
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B. High frequency dissipation: the JKD model 

The frontier between the low-frequency range (LF) and the high-frequency range (HF) is reached when the viscous 
efforts are similar to the inertial effects. The transition frequency is given bjii2 



fc 



r)(f> 



(5) 



2 7T a K pf 2 7T 

In LF, the flow in the pores is of the Poiseuille type, and dissipation efforts in the second equation of @ are given by 

F(t) = S(t) «=> F(t) * w(x, y, t) = w(x, y, t), (6) 

where d is the Dirac distribution. In HF, the width of the viscous boundary-layer is small in comparison with the size of 
the pores, and modeling the dissipation process is a more complex task. Here we adopt the widely-used model proposed 
by Johnson-Koplik-Dashen (JKD) in 1987, which is valid for random networks of pores with constant radiii^ The only 
additional parameter is the viscous characteristic length A. We take 



P : 



4 a k 
M2' 



P 



■qcj) 2 A 2 
4 a 2 k 2 pf 



(7) 



where P is the Pride number (typically P as 1/2). Based on the Fourier transform in time, F(td) = J R F(t)e lult dt, the 
frequency correction given by the JKD model can be written 



F{u) 



1 + i ui 



4 a 2 k 2 pf 
r\ A 2 <j) 2 

1/2 



1/2 



= I 1 + iP — 



= -=((!■ 



(8) 



a/2 



This correction is the simplest function satisfying the LF and HF limits of the dynamic permeability. 3 Therefore, the term 
F(t) * w(x, y, t) involved in the second equation of ([2]) is 



(9) 



The operator D 1 / 2 is a shifted order 1/2 time fractional derivative, generalizing the usual derivative characterized by 
= J 7 ' 1 (iuiw(x, y,u))). The notation (D + fi) 1/2 accounts for the shift fl in ©. 



C. The Biot-JKD equations of evolution 

Based on ([2]), f4]) and (j9]), the Biot-JKD equations can be written 



<9v s 



dt 
<9v s 



Pf 



dw 
~dt 



V<T, 



Pf 



dt 



^^ + ^^(Z? + f7)V 2 w = -V P , 

Ot K y/tt 



(10) 



a = (A/tre- j8m0l + 2jue, 
t p = m (— j3 tr £ + ^). 



We rearrange this system by separating and ^j- in the first two equations of (|10[) and by using the definitions of e 
and ^. Taking 



= V P 1 



(ii) 
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one obtains the following system of equations of evolution 
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(12) 



Terms f Vax , f Vsy , f Wx , f Wy , f axx , f„ xyi f ayy and f p have been to simulate sources. 



D. The diffusive representation 

Taking 



d w 

D Q w(x,y,t) = — + flw, (13) 
dt 



the shifted fractional derivative ((9]) can be written^* 

{D + Q) 1 / 2 w(x,y,t) = ^= f e D n w(x,y,r)dT. (14) 

V 71 " Jo Vt-T 

The operator [D + Q) 1 ' 2 is not local in time and involves the entire time history of w. As we will see in section UTTl 
a different way of writing this derivative is more convenient for numerical evaluation. Based on Euler's T function, the 
diffusive representation of the totally monotone function is^ 

1 1 f°° 1 

4= = -^/ "4 e ^- (15) 

Vt V* Jo Vo 

Substituting (JTSJ into <jT4j) gives 

1 f l f°° 1 

(D + SlfPwiwt) =~ ^ e ~ Ht ~ T) e- ^ D n w(x,y,T)dedT, 

7T Jo Jo VO 

(16) 

= - / —rP(x,y,0,t)de, 
n Jo VO 

where the diffusive variable is defined as 

i>(x,y,6,t)= [ e-WN^w^^jfr. (17) 
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For the sake of clarity, the dependence on and w is omitted in ■0. From (JTTJ) , it follows that the diffusive variable xp 
satisfies the ordinary differential equation 

**=-(* + n)V, + Al w, (i8) 

ip(x,y,0,0)=0. 

The diffusive representation therefore transforms a non-local problem (fT4|) into a continuum of local problems (fT8|). It 
should be emphasized at this point that no approximation have been made up to now. The computational advantages 
of the diffusive representation will be seen in sections Mil and [Vj where the discretization of ([TO)) and (fT5)) will yield a 
numerically tractable formulation. 



E. Energy of Biot-JKD 

Now, we express the energy of the Biot-JKD model (JTUJ). This result generalizes the analysis performed in the ID case 
in the Ref. [lfl 

Proposition 1 Setting C the 4x4 poroelastic tensor such that <j = C e — pi in 0), let 

E = Ei + E2 + E3 , 

with 



E 1 = \ f (pv 2 + p w w 2 + 2p f v s .w)dxdy, 



E 2 = 1 / [Ce : e + -p 2 ) dxdy, 

Iff V 1 1 I , 



Es = - / / J.--———(v- 1 l,yd0dxdy. 



Without any source terms, E is an energy which satisfies 
dE f f r] 1 1 I 



dt 



f f ll*-±— (nw 2 + (o + n)4, 2 )dedxdy < 0. (20) 

Jr2 Jem+ K * VtlO + 211 



Proposition [T] calls for the following comments: 

• the Biot-JKD model is well-posed; 

• when the viscosity of the saturating fluid is neglected (rj = 0), the energy of the system is conserved; 

• the terms in (fT9|) have a clearly physical significance: E\ is the kinetic energy, and E2 is the potential energy. The 
term E3 corresponds to the kinetic energy resulting from the relative motion of the fluid inside the pores. 



F. Dispersion analysis 



A plane wave de^ ajt_k r - 1 is injected in (fT2|) . where k = ke and d are the wavevector and the polarization, respectively; 
r is the position, uj — 2 7r / is the angular frequency and / is the frequency. If d is collinear with k, the dispersion relation 
of compressional waves is obtained. Setting 



f D A 



m (An + 2/z), 



D 2 M = - ((A/ + 2 fx) p w + m (p - 2 pf 0)) u? + iu-F(uj) (Xf + 2/i), 

K 

Dq{lu) = X u 4 -iuJ 3 - P F(uj), 

K 



(21) 
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the dispersion relation of compressional waves takes the form 

D e (k, uj) = D A k 4 + D 2 (uj) k 2 + D (uj) = 0. (22) 
If d is orthogonal with k, the dispersion relation of the shear wave is obtained. Setting 

' D 4 (uj) = lo 2 (p + cj) p f (a - 2)) - ioj (f 2 - F(u), 
D 2 (u>) =-oj 2 <t>pf{a-l) + ioj(f) 2 ^-F{uj), (23) 

K 

Dq(ui) — uj 2 <f) Pf a, — iui <j) 2 — F(u>), 

k 



the dispersion relation of shear wave takes the form 



1 {D A D Q -L><;\ 1/2 



k = — * " 2 . (24) 

Expressions (|2Tt - (|22p - ([2^)) - (IM)l are valid in the case of both the Biot-LF and Biot-JKD models with the frequency correc- 
tion defined by 

{F lf (uj) =1 Biot-LF, 
1 , ( 25 ) 

Fjkd(oj) =^=(0 + iw) 1/2 Biot-JKD. 

The solutions k p f, k ps of (|22p and the solution k s of (|24|) give the phase velocities c p / = w/5fte(fc p /) of the fast compressional 
wave, c ps = uj /$te(k ps ) of the slow compressional wave, and c s = w/5ie(fc s ) of the shear wave, with < c ps < c p f and 
< c s . The attenuations a p f = -3m(fc p j), a ps = — 3m(fc ps ) and a s = — 3m(fc s ) can also be deduced. Both the phase 
velocities and the attenuations of Biot-LF and Biot-JKD are strictly increasing functions of the frequency. The high 
frequency limits of phase velocities of compressional waves, c^j and c^, obtained by diagonalizing the left-hand side of 
system (|T2|) . satisfy the relation 



X c 4 - ((Xf + 2 ft) p w +m,(p-2p f l3)) c 2 +m (A + 2 p) = 0, (26) 
and the high frequency limit of phase velocity of the shear wave, cf , is 

c? = ^(p-^)~ 1 ' 2 . (27) 

Figure [1] shows the dispersion curves corresponding to the Biot-LF and Biot-JKD models. The physical parameters are 
those used in the numerical experiments presented in section fVl (medium VLq). Note that the scales are radically different 
in the case of fast and slow waves. The following properties can be observed: 

• when / < f c , the Biot-JKD and Biot-LF dispersion curves are very similar as might be expected, since 
lim F jkd (lj) = 1: 

uj— >0 

• the fast compressional wave and the shear wave are almost not affected by the frequency correction F(u>), while the 
slow compressional wave is greatly modified; 

• when / <gC f c , the slow compressional wave is almost static. When / > f C) the slow wave propagates but is greatly 
attenuated. 



III. THE BIOT-DA (DIFFUSIVE APPROXIMATION) MODEL 

The aim of this section is to approximate the Biot-JKD model, using a numerically tractable approach. 
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Figure 1. Dispersion curves: comparison between Biot-LF and Biot-JKD. 
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A. The Biot-DA first-order system 



Using a quadrature formula on N points, with weights ag and abscissa 9g > 0, the diffusive representation (|16[) can be 
approximated by 



f f°° 1 

{D + ny/ 2 w(x,y,t) = - / -=tl>(x,y,t,0)M, 
11 Jo V0 

N 

- 22aiip(x,y,t,6i), 

N 



(28) 



>From (|T8j) . the N diffusive variables if>g satisfy the ordinary differential equations 

= -(0 t + n,)il>t + Dnw, 
^e(x,y,0) = 0. 



(29) 



The fractional derivatives are replaced by their diffusive approximation (f2"5)) in the JKD model (|12p . Upon adding the 
equations (|29p and performing some straightforward operations, the Biot-DA system is written as a first-order system in 
time and space 



d v. 
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dp 



( da xx 
V dx 


! do X y 

dy 


( da xy 


f dCTyy 


\ dx 


dy 


f d(Txx 


da xy s 


\ dx 


dy t 


r da xy ^ 


dOyy' 


v dx 


dy / 
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p dp 
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p dp 



N 



^ w y ~ 7 ^ at 



ij -\- r^) ^p x j -\- fwxi 



dt x \ dx dy J x dy 
where j = 1 , N. Taking the vector of unknowns 

U = (v sx ,v sy ,w x , Wy, a xx , a X y,a yy ,p,ip x i,ipyi, . . . , ipxN, ^ v n) T 

and the source vector 

F = (/» IS 1 fv y3 1 fiUx 7 fwy 1 fo X x i fiTxy J j CT y y J ^ J f W ^ J fvjy ) • • • ) /u) x J /lOj, ) J 



(30) 



(31) 
(32) 
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the system ([50)) can be written 



<9U 9U dXJ 

7 - + A— +B— = -SU + F, 33 

at ox ay 

where A is the (2 N + 8) 2 propagation matrix along x, B is the (2JV + 8) 2 propagation matrix along y and S is the 
(2 7V + 8) 2 dissipation matrix l|A H IA2 ) lA3 ]l . The size of the system increases linearly with the number N of diffusive 
variables. 



B. Properties 

Four properties of system (|3"3")l are specified: 

• the eigenvalues of A (|A"Tj) and B (|A"2)l are real: with multiplicity 2N + 2, ±c™ f , ±c^ and ±c£°. The system 
is therefore hyperbolic; 

• since the eigenvalues and eigenvectors do not depend on the diffusive coefficients, they are the same in both the 
Biot-DA and Biot-LF or Biot-JKD models. This is not so in the case of the method presented in the Ref. 0, where 
the propagation matrix is modified to account for the fractional derivative; 

• the dispersion analysis is obtained in the case of the Biot-DA model by replacing F by 

! + l U! y—^ ag 

£rf Oe + n + iu 



o i ■ N 

F ^H = - 7 ^E^T- ( 34 ) 



in equations fl2"T|) -([2"2"l). 



C. Determining the Biot-DA parameters 

The coefficients an and 9t in ([2^)1 have to be determined. For this purpose, we want to approach the original frequency 
correction Fjkd(u>) (|25p by Fad{u) (I34[) in a given frequency range of interest. Let Q(to) be the optimized quantity and 
Q re f{to) be the desired quantity 



Fad(u) 
( W ) 



JKD 



N 

E 



(O + jw) 1 / 2 



ag ■ 

6t + + l Ld 

1 i=\ 



(35) 



k Qre/(w) = 1. 



We implement a linear optimization procedur o 10 ' 16 " — in order to minimize the distance between Q{uS) and Q re /(w) in the 
interval [uj m i n ,u> max ] centered on wo = 2tt fo, where /o is the central frequency of the source. The abscissas dg are fixed 
and distributed linearly on a logarithmic scale 



ig — U)„ 



!jj n 



e-i 

N-l 



£=1,...,N. 



(36) 



The weights ag are obtained by solving the system 

N 



}]atqt(u)k) = 1; k = l,...,K, 



(37) 



where the ujk are also distributed linearly on a logarithmic scale of K points 



uj k = u„ 



fc— i 

K — l 



k = 1 K. 



(38) 
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Since the qtioj) are complex functions, optimization is performed simultaneously on the real and imaginary parts 



/v 



^2a e Re(q e (uj k )) = 1, 
1=1 

N 

}^aeIm(qe(Q k )) = 0, k=l,...,K. 



(39) 



For practical purposes, we use uj m i n — wo/10 and ui max — IOujq, and we solve an overdetermined system (K = N), which 
can be solved by writing normal equations.— The number of diffusive variables N can be determined in terms of the ratio 
fo/fc and of the chosen accuracy: see Ref. [l0| for details. 



IV. NUMERICAL MODELING 
A. Splitting 

In order to integrate the Biot-DA system (|33p , a uniform grid is introduced, with mesh size A x = A y and time step A t . 
The approximation of the exact solution U(x, = i A x, y 3 = j A y, t n = n At) is denoted by U^- . If an unsplit integration 
of (|33[) is performed, a Van-Neumann analysis typically yields the stability condition 




A/ ,„n„[T— (40) 

where R(S) is the spectral radius of S, and T > depends on the numerical scheme. We have no theoretical estimate of 
R(S), but numerical studies have shown that this value is similar to that of the spectral radius in LF: 2 £ which can be 

very large£ The time step can therefore be highly penalized in this case (|40p. 

A more efficient strategy is adopted here, which consists in splitting the original system (|33ll into a propagative part 
and a diffusive part: 

— A— +B— -0 

dt dx dy ' 

au 

For the sake of simplicity, the source term F has been omitted here. The discrete operators associated with the propagative 
part and the diffusive part are denoted by H a and respectively. The second-order Strang splitting is then used to 
integrate |55]) between t n ant t n +i, giving the time-marching 



UW = H 6 (A*)U«, 

U\? = H a (At)ug } , (42) 



. U^ 1 = H 6 (A*) U g ) . 

The discrete operator H a associated with the propagative part (first equation of (|4ip) is an ADER 4 (Arbitrary 
DERivatives) schemed This scheme is fourth-order accurate in space and time, is dispersive of order 4 and dissipative of 
order 6, and has a stability limit T = 1. On Cartesian grids, ADER 4 amounts to a fourth-order Lax-Wendroff scheme. 

Since the physical parameters do not vary with time, the diffusive part (second equation of I41[) ) can be solved exactly. 
This gives 



At 
2 



K b [=r)U l3 =e-^ s U l3 . (43) 



The matrix e 2 s is computed numerically using the (6, 6) Pade approximation in the "scaling and squaring method"^ 
The numerical integration of the diffusive step (|43|) is unconditionally stabled 



The full algorithm (|4"2"]l is therefore stable under the optimum CFL (Courant-Friedrichs-Lewy) stability condition 

CFL = C ^ X"x < T ' (44) 
which is always independent of the Biot-DA model coefficients. 
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B. Immersed interface method 




1 



n 



Figure 2. Interface F between two poroelastic media Qo and Oi. 



Let us consider two homogeneous poroelastic media fio and fii separated by a stationary interface T, as shown in figure 
[H The governing equations in each medium have to be completed by a set of jump conditions. The simple case of 
perfect bonding and perfect hydraulic contact along T is considered here, modeled by the jump conditions: 22 



The discretization of the jump conditions requires special care, for three reasons. First, if the interfaces do not coincide 
with the uniform meshing, then geometrical errors will occur. Secondly, the jump conditions (|45[) will not be enforced 
numerically by the finite-difference scheme, and the numerical solution will therefore not converge towards the exact 
solution. Lastly, the smoothness of the solution required to solve (|33[) will not be satisfied across the interface, which will 
decrease the convergence rate of the ADER scheme. 

To overcome these drawbacks without detracting the efficiency of Cartesian grid methods, an immersed interface 
mcthod^ 2 ^ is used. The latter studies can be consulted for a detailed description of this method. The basic princi- 
ple is as follows: at the irregular nodes where the ADER scheme crosses an interface, modified values of the solution are 
used on the other side of the interface instead of the usual numerical values. 

Calculating these modified values is a complex task involving high-order derivation of jump conditions (|45[) . Beltrami- 
Michell relations ([3]) and singular value decompositions. Fortunately, all these time consuming procedures can be carried 
out during a preprocessing stage and only small matrix-vector multiplications need to be performed during the simulation. 
After optimizing the code, the extra CPU cost can be practically negligible, i.e. lower than 1% of that required by the 
time-marching procedure. 

The immersed interface method uses only the propagative part (first equation of (|4"Tj) ) and the jump conditions (|4"5")) . 
whose do not depend on the frequency regime. Consequently, no changes are required in the high-frequency range compared 
to the low-frequency range. The algorithm is the same as in the low-frequency range, see Ref. [H for details. 

V. NUMERICAL EXPERIMENTS 

General configuration 

The physical parameters used in all the numerical experiments, which are given in table HI correspond to Cold Lake 
sandstone (CIq) and Berea sandstone (Qi) saturated with water. The computational domain [— 0.1,0. 1] 2 m is discretized 
with N x x N y grid points, N x — N y = 1600, which amounts to 29 points per slow compressional wavelength in f2o an d in 
fii for a source of central frequency fa = 200 kHz. As seen in ID, 10 N = 6 diffusive variables are used to ensure an error 
of model smaller than 5.58 %. The time step is deduced from (|4"4")l . taking CFL = 0.95. The numerical experiments are 
performed on an Intel Core i7 processor at 2.80 GHz. 



[v.] = 0, [w.n] = 0, [<r.n] = 0, [p] = 0. 



(45) 



Biot-JKD poroelastic waves in 2D heterogeneous media 11 





Parameters 


Ho 




Saturating 


; fluid p f (kg/m 3 ) 


1040 


1000 




n (Pa.s) 


1.5 10~ 3 


10" 3 


Grain 


p s (kg/m 3 ) 


2650 


2644 




p (Pa) 


2.93 10 9 


7.04 10 9 


Matrix 




0.335 


0.2 




a 


2 


2.4 




k (m 2 ) 


io- n 


3.6 10~ 13 




\ f (Pa) 


6.14 10 9 


1.06 10 10 




m (Pa) 


6.49 10 9 


9.70 10 9 




P 


n Kt? in — 1 
9.00 10 


!\2U 10 




A(m) 


2.19 10~ 5 


5.88 10" 6 


Phase velocities djff (m/s) 


2384.17 


3269.89 




(m/s) 


758.95 


814.95 




c? (m/s) 


1229.00 


1776.16 




/c (Hz) 


3.84 10 3 


3.68 10 4 



Table I. Physical parameters used in numerical experiments. 



Test 1: homogeneous medium 

The unbounded homogeneous medium fin is excited by a point source. The only non-null component of the vector of 
source F (j3"2")) is / CTxb = g(x, y) h(i), where g(x, y) is a truncated gaussian centered at point (0, 0), of radius i?n = 3.79 10~ 3 
m and E = 1.90 HT 3 m: 

\ ( x 2 + v 2 \ 

exp -JL X0^x 2 +y 2 ^R 2 , 




g(x,y) = { ttE* ' V »' 7 " i4(» 

otherwise, 

and ft(t) is a Ricker signal of central frequency /q = 200 kHz and of time-shift to = 2//o = 10~ 5 s: 

if0<t<t , (47) 
otherwise. 

We use a truncated gaussian for (7(2, y) rather than a Dirac to avoid spurious numerical artifacts localized around the 
source point. This source generates cylindrical waves of all types: fast and slow compressional waves and shear wave, 
which are denoted Pf, P s and 5, respectively, in figure [3] Fast and slow compressional waves are observed as regards the 
pressure, while the additional shear wave is present in the a yv component of the stress tensor. No special care is applied 
to simulate outgoing waves (with PML, for instance), since the simulation is stopped before the waves have reached the 
edges of the computational domain. 

Test 2: diffraction of a plane wave by a plane interface 

In all the following tests, the source is a plane right-going fast compressional wave, whose wavevector k makes an angle 
6 = with the horizontal a:-axis. Its time evolution is the same Ricker signal as in the first test (|47|) . We use periodic 
boundary conditions at the top and at the bottom of the domain. 

This test is conducted on a vertical plane interface at x — m. The incident Pf -wave (Ip/) propagates in the medium 
£7o- The incident wave crosses the interface normally, leading to a ID configuration. Consequently, only the fast and 
slow compressional waves propagate. From a numerical point of view, however, the problem is fully bidimensional. The 
advantage of such a ID configuration is that each diffracted wave has interacted with the interface and is consequently 
very sensitive to the discretization of the interface conditions P5|) . The figure U shows a snapshot of the pressure 
at t\ = 1.82 10 -5 s, on the whole computational domain. The reflected fast and slow compressional waves, denoted 
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respectively Rpy and Rp s , propagate in the medium f2 ; and the transmitted fast and slow compressional waves, denoted 
respectively Tp/ and Tp s , propagate in the medium Cl±. 



(a) 



(b) 




-0.05 -0.04 -0.02 0.02 0.04 O.0C 0.03 

x(m) 



Tps I 



Tpf 



-0.05 -0.04 -0.02 0.02 0.04 O.0C 0.03 

x(m) 



Figure 4. Test 2. Snapshot of p at initial time (a) and at ti = 1.82 10 5 s (b). 



In this case, we compute the exact solution of Biot-DA thanks to standard tools of Fourier analysis. Special care is 
taken to ensure that the initial data and the exact solution are valuable reference solutions : Nf = 512 Fourier modes are 
used, with a frequency step A / = 4000 Hz. The figure [S] shows the excellent agreement between the analytical and the 
numerical values of the pressure along the line y = m. 



Test 3: diffraction of a plane wave by a cylinder or a circular shell 

In the first case, a cylindrical scatterer Hi of radius 0.015 m is centered at point (0.01,0). In the second case, we 
consider a poroelastic shell Sli of external radius 0.03 m and internal radius 0.015 m, centered at point (0.025,0). The 
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Pressure 



Zoom on the slow compressional waves 



800 



600 



400 



200 

Q_ 



-200 



-400 
-0 



analytical 

o numerical 



-0.05 





0.02 -0.015 -0.01 -0.005 0.005 0.01 0.015 0.02 
x (m) 



Figure 5. Test 2. Pressure along the line y — m; vertical line denote the interface. Comparison between the numerical values 
(circle) and the analytical values (solid line) of p at t\ = 1.82 10~ 5 s. 



initial conditions are illustrated in figure [6fa) and Ela) , and the figure |6][b) and [TJb) represent a snapshot of a yy at 
<i = 3.63 10 -5 s. The immersed interface method ensures that no spurious diffractions are created during the interaction 
of the incident wave with the scatterers, despite the non-null curvature of the interfaces. In the figure [6jb) and[7][b), the 
transmitted fast compressional wave has a curved wavefront because its phase velocity c p / in the medium Ox is greater 
than in the medium f2o- Classical conversions and scattering phenomena are observed. 



(a) 



CM 
Ci, qs - 
0.05 " 
0.04 



-0.03 -0.05 -0.04 -0.02 0.02 0.04 0.05 0.03 0.1 

x(m) 




-O.0S -0.05 -0.04 -0.02 0.02 0.04 0.05 0.03 0.1 

x (m) 



Figure 6. Test 3. Scatterer: cylinder. Snapshot of a yy at initial time (a) and at ti = 3.63 10 5 s (b). 



Test 4-' multiple ellipsoidal scatterers 

In the last test, the ability of the proposed numerical strategy to handle complex geometries with variable curvature is 
illustrated. 15 ellipsoidal scatterers, of major radius 0.008 m and of minor radius 0.005 m, of medium f2j are randomly 
embedded in Qq. Performing such a simulation with hundreds of scatterers has physical applications. Doing so provide a 
mean to check the validity of multiple scattering models i 25 ' 26 
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VI. CONCLUSION 

An explicit finite-difference method is presented here for simulating transient poroelastic waves in the full range of 
validity of poroelasticity. The Biot-JKD model, which involves order 1/2 fractional derivatives, was replaced here by 
an approximate Biot-DA model, which is much more tractable numerically. The Biot-DA coefficients are determined 
here using an optimization procedure, which depends on the frequency range of interest. The hyperbolic system of partial 
differential equations was discretized using efficient tools (Strang splitting, fourth-order ADER scheme, immersed interface 
method). This numerical methodology provides highly accurate numerical solutions, allowing a fine numerical investigation 
of realistic porous media. 
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Some future lines of research are suggested: 

• Continuously variable porous media. In the previous numerical experiments, the physical parameter characterizing 
the porous media are constant and discontinuous across the interfaces. The possibility of applying the presented 
numerical method to 2D porous media with continuously variable coefficients, r» where no analytical expressions are 
available, is currently in progress: see Ref. [l0| in ID. 

• Multiple scattering. The numerical tools presented here make possible the modeling of multiple scattering in random 
media. Based on simulated data, the properties of the effective medium amount to the disordered under study can 
be deduced] 25 ' 26 Optimization and parallelization of the method must be realized to take into account of hundreds 
of inclusions in our simulations. 

• Anisotropy. The anisotropy of some porous media is sometimes a main feature, for instance in biomechanics to 
model trabecular and cortical bonesj 24 ' 27 The anisotropy of a medium changes only the 4x4 poroelastic tensor C 
([1]). So the ADER 4 scheme and the immersed interface method have to be modified, but the approximation of the 
fractionnal derivative is still valid. 

• Thermic boundary-layer. In cases where the saturating fluid is a gas, thermo-mechanical effects have to be taken 
into account. Extended versions of the Biot-JKD have been developed, 2 ® involving additional order 1/2 fractional 
derivatives. The numerical method developed in this paper should lend itself well to working with this model. 



Appendix A: MATRICES OF PROPAGATION AND DISSIPATION 

The matrices in 
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